Thermodynamic and Structural Properties of the High Density Gaussian Core Model 
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We numerically study thermodynamic and structural properties of the one-component Gaussian 
core model (GCM) at very high densities. The solid-fluid phase boundary is carefully determined. 
We find that the density dependence of both the freezing and melting temperatures obey the asymp- 
totic relation, log!/, logT m oc — p 2 / 3 , where p is the number density, which is consistent with Still- 
inger's conjecture. Thermodynamic quantities such as the energy and pressure and the structural 
functions such as the static structure factor are also investigated in the fluid phase for a wide range 
of temperature above the phase boundary. We compare the numerical results with the prediction of 
the liquid theory with the random phase approximation (RPA). At high temperatures, the results 
are in almost perfect agreement with RPA for a wide range of density, as it has been already shown 
in the previous studies. In the low temperature regime close to the phase boundary line, although 
RPA fails to describe the structure factors and the radial distribution functions at the length scales 
of the interparticle distance, it successfully predicts their behaviors at shorter length scales. RPA 
also predicts thermodynamic quantities such as the energy, pressure, and the temperature at which 
the thermal expansion coefficient becomes negative, almost perfectly. Striking ability of RPA to pre- 
dict thermodynamic quantities even at high densities and low temperatures is understood in terms 
of the decoupling of the length scales which dictate thermodynamic quantities from the interparticle 
distance which dominates the peak structures of the static structure factor due to the softness of 
the Gaussian core potential. 



PACS numbers: 



I. INTRODUCTION 



Complex fluids such as colloidal suspensions and emul- 
sions are often regarded as macroscopic models of atomic 
or molecular systems. They are ideal benches to test liq- 
uid theories developed to describe thermodynamic, dy- 
namic, and structural properties of atomic and molecular 
liquids jT|. It is not only because the size of constituent 
unit of complex fluids arc much larger than atomic coun- 
terpart but also because their interparticle interactions 
can be tailor-made and tuned relatively easily. While 
the pair interactions of atomic systems are exclusively 
characterized by short-ranged and strong repulsions with 
weak and longer-ranged attractions, leading to the typ- 
ical phase diagram demarcating gas, liquid, and crys- 
talline phases [l(, those for the complex fluids are far 
more diverse. This diversity leads to very rich and of- 
ten counter- intuitive macroscopic behaviors Q • Amongst 
those interactions, the ultra-soft potential have attracted 
particular attention recently in the soft- condensed mat- 
ter community [H-Hll- The ultra-soft potentials are 
isotropic repulsive potential characterized by weak and 
bounded repulsion at short distance and the mild re- 
pulsive tails whose steepness is much smaller than the 
typical atomic potential. These potentials are realized 
in many complex fluids such as star-polymers (HHHI, 
dendrimers H, H(|, and the polymers in good sol- 
vent @, l3lM34 1 . Thermodynamic phase diagrams of the 
ultra-soft particle systems have very distinct and exotic 
properties such as the re-melting from solid to fluid phase 
at high densities, the re-entrant peak at the intermedi- 
ate densities @, i, i, GH, HH, HI HI, negative thermal 



expansion coefficient [H, [T(| , and the cascades of the var- 
ious crystalline phases at very high densities J^HH- The 
Gaussian core model (GCM) is one of the simplest exam- 
ples of the ultra-soft potential systems. GCM consists of 
the point particles interacting with a Gaussian shaped 
repulsive potential; 



v{r) = eexp[— {rja)' 



(1) 



where r is the interparticle separation, e and a are the pa- 
rameters which characterize the energy and length scales, 
respectively. GCM was first introduced by Stillinger Q 
and has been studied by many groups [5Hla|. Despite of 
its simple form of the pair potential, GCM exhibits many 
typical thermodynamic behaviors of the ultra-soft par- 
ticles. According to thermodynamic phase diagram ob- 
tained by numerical simulations [3, HJ , GCM basically be- 
haves like hard spheres at very low densities and temper- 
atures; the crystalline structure in the solid phase is fee 
and the freezing/ melting temperatures sharply increase 
with density. However, as the density increases further, 
the freezing temperature reaches a maximal value and 
beyond this point it changes to a decreasing function of 
the density. This re-entrance takes place at per 3 0.25, 
where p is the number density. Concomitantly, the crys- 
talline structure changes from fee to bee. Recently, ther- 
modynamic and transport anomalies of the fluid phase 
in the vicinity of the reentrant peaks are investigated 
[Tol . [l2rfl6f • Microscopic and structural properties such 
as the static structure factor in the fluid phase are also 
reported and documented [1, 0, HI, UK ■ These studies re- 
vealed that, as the density increases beyond the reentrant 
peak but at a fixed temperature, thermodynamic and 
structural properties of GCM becomes more ideal-gas- 
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like, signaled by the lowering of the peak of the structure 
factors and the better agreement with simple approxima- 
tions such as the random phase approximation (RPA). 
Most of studies in the past, however, have focused on the 
densities not far from the reentrant peak or the relatively 
high temperatures. Less attention has been paid for the 
high density and low temperature regimes, especially in 
the vicinity of the solid-fluid phase boundary. Near the 
phase boundary line, the thermodynamic and structural 
properties are expected to be highly non-trivial even at 
the high density limit. Based on the duality argument of 
the ground state of GCM in the reciprocal space, Still- 
inger has conjectured that the freezing and melting tem- 
peratures, Tf and T m , are given by an asymptotic form; 

logT/, logT m cx V /3 (2) 

in the high density limit However, this conjecture has 
not been confirmed numerically. Recently, we studied the 
nucleation and glassy dynamics of the one-component 
GCM in the supercooled state at the unpreccdcntcdly 
high densities, 0.5 < pa 3 < 2 [35[ . It was found that the 
crystal nucleation rate decreases drastically as the den- 
sity increases and concomitantly dynamics of the con- 
stituent particles becomes very sluggish. The density 
time correlation function exhibits typical behavior of the 
supercooled liquids near the glass transition point, such 
as the two-step and non-exponential structural relax- 
ation. The relaxation time steeply increases as the tem- 
perature is lowered at a fixed density. Surprisingly, these 
are well described by the mode-coupling theory, implying 
that the high density and one-component GCM is more 
amenable to the mean-field picture of the glass transition 
than other typical glass formers. These observations call 
for more detailed analysis of the high density GCM at 
the low temperature regime. Especially, it is tempting 
to consider GCM in the high density limit as the ideal 
and clean model system to study the glass transition. 
Thermodynamic and structural characterization are pre- 
requisites for dynamical study [HI, [3(| but the detailed 
study is lacking. 

In this work, we numerically investigate thermody- 
namic and structural properties of the one-component 
GCM up to the density pa 3 = 2.4. We determine the 
solid-fluid phase boundary and show that the Stillinger's 
scaling, Eq. ©, holds at pa 3 > 1.2. Thermodynamic 
and microscopic structural properties of GCM are also 
analyzed carefully over a wide range of temperature and 
density. The potential energy, pressure, thermal expan- 
sion coefficient, and the static structure factors are evalu- 
ated and compared with the prediction of the liquid state 
theory. Surprisingly good agreement with the random 
phase approximation (RPA) is found for thermodynamic 
quantities for a wide range of temperature, including the 
low temperature regimes where the same approximation 
poorly describes the static structure factor and radial 
distribution function. This counterintuitive observation 
can be attributed to the ultra-soft nature of GCM for 
which the microscopic structure near the first shell of the 



system decouples with the macroscopic properties. 

This paper is organized as follows. In Section II, tech- 
nical details of simulations and the method to compute 
the phase boundary are discussed. In Section III, we 
present simulation results for the phase diagram, thermo- 
dynamic quantities, and structural functions. We com- 
pare the simulation results in the fluid phase with RPA 
predictions in Section IV. We summarize and discuss the 
results in Sec. V. 



II. SIMULATION METHOD 

A. MD and MC simulation 

Thermodynamic state of GCM is fully characterized by 
the density p and temperature T. In this work, we focus 
on the density and temperature range of 0.3 < p* < 
2.4 and 10~ 6 < T* < 1, where p* = pa 3 and T* = 
k^T /e. In order to analyze thermodynamic properties 
and determine the solid-fluid phase boundary, molecular 
dynamics (MD) simulations with Nose thermostat are 
carried out under the periodic boundary conditions. To 
integrate the equations of motion, we use a reversible 
algorithm similar to the Velocity- Verlet method [32| with 
time steps of O.lr, which is sufficientl y short to preserve 
the Nose Hamiltonian. Here r = ^Jma 2 Je is the time 
unit, where m is the mass of a particle. For evaluation of 
the free energy of the reference state (see below), Monte 
Carlo (MC) simulation is used. In a trial MC move, the 
maximum displacement of a particle is adjusted to keep 
the acceptance ratio about 50 %. In both simulations, 
the total number of particles is N = 3456. This is twice 
the cube of an integer (in this case 12), a natural choice 
for the bec crystal in a cubic simulation box. The cutoff 
length of the potential is taken as 5cr. The pair potential 
at the cutoff length is 1.4 x 10 _11 e which is much smaller 
than the typical kinetic energy at the lowest temperature 
studied in this work. 



B. Evaluation of the free energy 

The chemical potential as a function of the tempera- 
ture and pressure, n{T, P), is required in order to deter- 
mine the phase boundary. We evaluate it using the free 
energy f(T 7 p) and pressure P(T, /?). 

We calculate the free energy using the thermodynamic 
integration method combined with the particle inser- 
tion method I38| for the fluid phase and the Frenkel- 
Ladd method [39J for the crystalline phase. This pro- 
cedure is the same as the one employed by Prestipino et 
al. H for lower densities. The free energy of the sys- 
tem is the sum of the ideal fid and excess part f ex . 
fid is given by f ld {T 1 p) = k B T{logA 3 pa 3 - 1), where 
A = \J 2-kH 2 / rak^T is the de Broglie thermal wave length. 
According to the thermodynamic integration scheme, 
f ex (T, p) can be evaluated by integrating over the energy 
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and pressure from the reference state point (To, po) to the 
target state point (T, p) using the following equations. 



fex(T, po) _ fex(T ,p ) 

T T 

fex(T, p) = fex(T, Po) 
T T 



T dT> < T '^ 

To 



JV2 



P dp/ fP(T,P') 1 



Po 



p' 2 T 



P' 



(3) 



where u is the potential energy per particle. For the fluid 
phase, the reference free energy fex(To, pp) is calculated 
by the particle insertion method [371 l38l ] and the pressure 
is evaluated from the virial equation. In this method, the 
free energy is expressed in terms of the energy cost to 
insert one particle into the system as 



fex(To,Po) 

To 



-fc B log ycxp 



knTo 



Po 
PoTa 



1,(4) 



where Ei n is the interaction energy of an inserted particle 
with other particles in the system. The average should be 
taken over the ensemble of randomly inserted particles. 
For the crystalline phase, on the other hand, the ref- 
erence free energy is computed using the Frenkel-Ladd 
method [13, Hi| which is a different kind of thermody- 
namic integration technique. In this method, we consider 
a hybrid Hamiltonian V(X), which interpolates between 
the Hamiltonian of the original system V and that of the 
Einstein crystal V e in as V(X) = V + (1 — X)V e i n , where A 
is the switching parameter. The free energy of the orig- 
inal system can be computed by the following equation, 
which is the integral over A of the Hamiltonian of the 
hybrid system evaluated from the simulation, 



J ex fe . 



v ( 



em I \ ! 



(5) 



where is the ensemble average under a hybrid 

Hamiltonian V^(A) and f ex ,ein is the excess part of the 
free energy of the Einstein crystal. We choose (T *,Pq) 
=(0.1, 0.01) and (T*, p*) = (0.0794, 0.28) as the ref- 
erence states for the fluid and crystalline phases, respec- 
tively. The ensemble averages in Eqs. (j4|) and ([5]) are eval- 
uated using the MC simulations. The integration over A 
in Eq. (J5J) is calculated by slicing A to the grids of the 
width of 0.05 and evaluating the value of the integrand at 
each grid point from independent MC simulations. Like- 
wise, the integral over the isothermal and isochoric path- 
ways in Eq. @ is computed by slicing the pathways into 
many grid points. The energy and pressure at each grid 
point are computed using the MD simulations. The free 
energy for both the fluid and crystalline phases are ob- 
tained by combining these data points and the reference 
free energy. In order to determine the solid-fluid phase 
boundary over the density range of 0.3 < p* < 2.4 with 
satisfactory accuracies, more than 800 grid points were 
necessary. 
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FIG. 1: The freezing (T?) and melting (T^) temperatures of 
GCM as a function of density (open up/down triangles). The 
result of Prestipino et al. is also plotted (filled circles with 
short-dashed line) The long-dashed line is the threshold 
temperatures, Tth, above which RPA gives a reasonable de- 
scription of the system (see Sec. IV). Open squares and dot- 
ted line indicate the temperature below which the thermal 
expansion coefficient becomes negative, T a , obtained from 
simulation and RPA, respectively. The short-dashed line at 
p* ~ 0.15 demarcates the fee (left) and bee (right) crystalline 
phases. 



III. SIMULATION RESULTS 

A. Phase diagram 

In his pioneering work, Stillinger has conjectured that 
the solid-fluid phase boundary is asymptotically given by 
Eq. ([2]) in the high density limit (4J. This conjecture is 
based on analysis of the density dependence of the po- 
tential energy of various crystalline structures (bec, fee, 
etc) at T = 0. According to his analysis, the potential 
energy is expressed as 

u=- e - + ^f^{l + AcM-Kp^) + ■■■}, (6) 

where A and K are constants which depend on the crys- 
talline structure. Stillinger argued that this density de- 
pendence remains qualitatively unchanged at finite tem- 
peratures. Since the ordered structure of the crystalline 
phase is responsible for the term exp(— Kp 2 ' 3 ) in Eq. ([6]), 
the energy difference between the crystalline and fluid 
phase at the phase boundary should be also proportional 
to this term. This argument leads us to a conjecture that 
the melting/freezing temperature is proportional to this 
factor, which is Eq. Here we verify this argument 
numerically. 

Figure [1] presents the phase diagram of GCM obtained 
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FIG. 2: -logT„ versus p* 2/3 . The solid straight line is a 
guide for the eyes. 



from our simulation. Both the freezing and melting tem- 
perature Tj? and are shown in this figure, but their 
values are very close to each other and indistinguishable 
in the scale of the figure. As expected, the melting and 
freezing temperatures dramatically decrease as the den- 
sity increases, down to T* w at the highest density 
we studied. The phase boundary for p* < 0.7 obtained 
by Prcstipino et al. Q is also plotted, in order to con- 
firm that the present result perfectly matches with theirs 
for the density window where both results arc available. 
The crystalline structure at high densities is bec, as has 
been verified by the direct MD simulation of the nucle- 
ation [35l . l36j . In order to verify the scaling relation, 
Eq. ([2]). we plot the logarithm of the melting tempera- 
ture as a function of p 2 / 3 in Fig. [2] One observes that 
the result rides on the scaling function at p* > 1.2. 



B. Potential energy and pressure 
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FIG. 3: The potential energy difference between the crys- 
talline and fluid phases against the melting temperatures with 
the fit by a straight line of the slope I (solid line). 



The first order transition from the crystalline to fluid 
phase is accompanied with the discontinuous change 
in the structural order and the potential energy. We 
show the temperature dependence of the potential en- 
ergy difference between the two phases u*f luid (Tf , p) — 
u* crystal {Tf,p) as a function of T m in Fig. where 
u* = u/e is the dimensionless potential energy. This 
figure shows that the energy difference is proportional to 
the melting temperature at the low temperatures/high 
densities, verifying the assumption which Stillinger has 
employed to conjecture Eq. The entropy differ- 

ence between the two phases can be estimated from 
this energy difference by s fluid (T f ,p) - s crystal (T f , p) = 
{ufi u id{Tf,p)—Ucrystal(Tf,p)}/Tf. From the result 
of Fig. [31 the entropy difference at the low tempera- 
ture/high density limit can be estimated as 

sjiuid(Tf,p) - s cryst ai(Tf , p) ~ 0.45fc B . (7) 

This value should be compared with the results at lower 
densities in the earlier work; 0.81fce at p* = 0.4 and 
0.54fc B at p* = 1.0 0. 
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FIG. 4: The equation of state of GCM. Filled and open 
circles are the results for fluid and crystalline (bec) phases, 
respectively. Solid lines are the results from RPA (see Sec. 
IV). Left panels are P-p plots at T* = 7.94 x fCT 3 (a), 5.0 x 
fCT 4 (c), and 7.9 x 1CT 6 (e). The inset in (a) is the closeup 
around the freezing transition density. Right panels are P-T 
plots at p* = 0.33 (b), 0.99 (d), and 2.0f (f). and T* (see 
text) are indicated by arrows. 

Next we focus on the equation of state (EOS) of the 
high density GCM, i.e., the pressure as a function of 
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the density and temperature. The dimcnsionless pres- 
sure P* = Pa 3 /e is plotted in Figure 2J In Figures 2] 
(a), (c), and (e), we plot the isothermal cut of EOS 
and the isochoric cut in Figures 0] (b), (d), and (f). 
The parameters T* and p* for each adjacent figures 
have been chosen so as for them to share the common 
freezing points; Fig. @] (a) and (b) share the freezing 
point (Tt,p* f ) = (7.94 x 1CT 3 ,0.33), (c) and (d) share 
(T*,p* f ) = (5.0 x 1(T 4 ,1.00), and (e) and (f) share 
(Tj;,p* f ) = (7.9 x 1(T 6 ,2.01). Figure H (a) shows that 
the melting of the bcc crystalline phase (white circles) to 
the fluid phase (filled circles) takes place at at p* « 0.33. 
The inset shows the narrow coexistence region around the 
transition density, at which the pressure becomes con- 
stant. Similar behaviors of the first order transition are 
also observed in Figs. @] (c) and (c) but at much higher 
densities. 

Figure [4] (b) shows the equation of state at p* = 0.33 
over the temperature range of 1.0 x 10~ 3 < T* < 1.0. At 
this relatively low density, the pressure is a monotonically 
increasing function of the temperature, which is usual be- 
havior of ordinary fluids. However, at higher densities, 
as shown in Fig.0](d) and (f), there exists the tempera- 
ture regime in which the pressure becomes a decreasing 
function of temperature. In this regime, the thermal ex- 
pansion coefficient a — V~ 1 (dV/dT)p becomes negative. 
We determine the threshold temperature T a at which a 
changes its sign by fitting the pressure by a smooth poly- 
nomial function and plot T a in Fig. Q] (open squares) . At 
low densities, T a is located in the vicinity of the phase 
boundary. With increasing density, however, the differ- 
ence between T a and T m increases monotonically. Al- 
though the existence of the anomalous negative thermal 
expansion coefficient of GCM has been reported in the 
literatures [3, H| , its high density behavior has not been 
explored. We shall discuss this result and its asymptotic 
behavior of T a at high densities in the following section. 



CO 




FIG. 5: The static structure factors in the fluid phase at the 
freezing temperatures for p* = 0.33 (O); 1-00 (□), 1.65 (A), 
and 2.01 (v)- The dotted lines are the results obtained by 
the Fourier transformation of g(r). The solid lines are the 
results of RPA (see Sec. IV) . The inset is the semi-log plot of 
the main figure for S(k) < 10 _1 at low fc's. 



IV. RANDOM PHASE APPROXIMATION 
ANALYSIS 

A. Random phase approximation 

It is known that the thermodynamics and microscopic 
structure of the GCM fluid in high densities and temper- 
atures are described by the random phase approximation 
(RPA) remarkably well [H, Q . RPA is one of the approx- 
imation schemes of the liquid theory and a kind of the 
mean-field theory for the thermodynamics and structure 
of liquids [l[ . In this section, we discuss how this approx- 
imation works at much lower temperatures. 

It is convenient to divide thermodynamic quantities 
into uniform and fluctuation parts as follows. The po- 
tential energy can be represented as 



Static structure factor 



ir 3 / 2 pa 3 e 



(8) 



In Figure O we show the static structure factors S(k) 
in the fluid phase just above the freezing temperatures. 
S(k) obtained from the Fourier transformation of the 
radial distribution function g(r) is also shown. At all 
densities, S(kys exhibit sharp peaks similar to those of 
the ordinary simple fluids, such as the hard sphere and 
Lennard- Jones fluid, near the freezing temperatures. The 
peak position fc max shifts to high fc's as the density in- 
creases, fc max oc p 1 / 3 , as one should expects at high den- 
sities. Note that, however, the height of the first peak is 
about 3.1 at all densities, which is slightly higher than 
the universal value 2.85 for the ordinary fluids which is 
known as Hansen- Verlet criterion [l|, |4(J . 



where the first two terms are the uniform part and Au is 
the fluctuation part. Am can be expressed as 



where 



Au 



v{k) 



1 



k 2 i{k)S{k) dk, 



3/2 



e<r 3 exp(-A;V/4) 



(9) 



(10) 



is the reciprocal expression of v(r) . Likewise, the pressure 
can be written, using the virial equation, as 



P = k B Tp + 



ir 3 / 2 p 2 a 3 e 



AR 



(11) 
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where the first two terms are the uniform part and the 
third is the fluctuation part which can be written as 



AP 



9 
4tt 2 



(fc 2 - k i a 2 /6)v{k)S(k) dk. (12) 



In RPA, the direct correlation function of the system is 
approximated as crpaM = — /3u(r), where = l/k^T. 
This approximation makes it possible to express various 
static quantities in simple and analytic forms. The static 
structure factor can be expressed as 



•S'RPA(fc) = 



1 



l + pj3v(k)' 



(13) 



The fluctuation parts of the potential energy and pressure 
are expressed as Q: 



Aurpa 



■2l 

PC 



Li 3 /2(-7): 



AFrpa = -— |Li 3 / 2 (-7) 



Li5/2(-7)} > 



(14) 



where 7 = 7r 3 / 2 po^e/k^T is a dimensionlcss coupling 
parameter and Li t/ (x) is the ^-th polylogarithm func- 
tion 0,0- Furthermore, the radial distribution function 
at r = can be expressed analytically as 



ffRPA^ = 0) = 1 + 



1 



1-3/2 



per 



3 Li 3/2( 



- 7 ). 



(15) 



The second term on the right hand side of this expression 
is negative for arbitrary densities and temperatures. At a 
very low temperature, the modulus of this term becomes 
larger than the first, leading to an unphysical negative 
g(r = 0). Wc refer to this temperature as the threshold 
temperature T t h- Wc plot T t h in Fig.Q](long-dashcd line). 
At high densities, Tth can be expressed analytically as 



knT t 



th 3/2 3 

— = ir ' pa exp 



2/3 



(16) 



which is obtained by the asymptotic expansion of poly- 
logarithm function of Eq. (jT5j) (sec Appendix). Interest- 
ingly, the threshold temperature follows the same asymp- 
totic scaling law log T t h oc — p 2 / 3 as the melting and freez- 
ing temperatures, Eq. ([2]). Note that this asymptotic 
expression is very accurate down to moderate densities 




FIG. 6: The radial distribution function in the high tem- 
perature regime for (a) p* = 0.10, (b) 0.33, (c) 1.00, and 
(d) 3.00. Circles and lines are simulation and RPA results, 
respectively. Two results in each panel correspond to the 
results at the threshold temperature (T* = 0.823, 0.535, 
0.159, and 0.689 x 10 -2 ) and higher temperature at which 
9RPA(r = 0) = 0.83 (T* = 5.76, 5.35, 4.30, and 2.29). 



We also computed the radial distribution function g(r). 
Fig.|5]shows g(r) obtained from simulation (open circles) 
and RPA (solid lines) at T = T t h and at much higher 
temperatures at which <7rpa(?" = 0) = 0.83, for several 
densities. At high temperatures, agreement of simula- 
tion results with RPA is excellent in all densities and 
for all r's. At T = Tth, however, RPA works poorly 
around r = 0, as expected from Eq. (fT5|) . On the other 
hand, RPA's results perfectly match with the simulation 
results at larger r's including the first shell peaks. The 
reason why thermodynamic quantities are well described 
by RPA even at Tth whereas agreement of g(r) near r = 
is poor can be attributed to the fact that the short-range 
part of g(r) does not contribute to both the potential 
energy and pressure, as wc shall discuss in the following 
subsection. 



B. High temperature regime 

We first assess the validity of RPA at temperatures 
above Tth- Tabic U shows the ratio of RPA to simulation 
results of the fluctuation parts of the potential energy and 
pressure at T = T t h- The deviations are smaller than 10% 
and they monotonically become smaller as the temper- 
ature increases, implying the thermodynamic quantities 
are well-described by RPA even at T t h- 



TABLE I: RPA results compared with simulation results for 
the fluctuation parts of the potential energy and pressure at 
the threshold temperatures for various densities. 



p 


0.10 


0.33 


1.00 


3.00 


Tth 


0.823 


0.535 


0.159 


0.689 


Amrpa/Au 


0.97 


0.94 


0.94 


0.98 


APrpa/AP 


1.07 


0.96 


0.92 


0.94 
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C. Low temperature regime 

We move to temperatures below T t h and discuss the 
validity of RPA in describing thermodynamic properties 
of GCM at high densities. In Fig. [51 the static struc- 
ture factors obtained from RPA, S , rpa(/c), are shown in 
the solid lines. It is obvious that RPA cannot capture 
even the qualitative behaviors of S(k); SrpaCO remains 
flat at higher fc's and does not possess any prominent 
peak. However, as the main panel and inset of Fig. [5j 
shows, RPA correctly predicts the low k's behavior up to 
just below the wavevector at which the first peaks are lo- 
cated. It is in stark contrast with ordinary simple atomic 
fluids for which RPA works poorly for the whole range of 
wavevectors [43j]. The excellent agreement implies that 
the mean-field character of the dense GCM still survives 
in the length scales slightly longer than the typical inter- 
particle distance. 

Next, we compare RPA with simulation results for 
thermodynamic quantities. We first look at the equa- 
tion of state. In Fig. [H the pressure obtained from RPA 
are plotted in solid lines. As shown in Fig. [4] (a) and (b), 
the deviation of the values of RPA from those of simu- 
lation is very large at p* = 0.33 and the discrepancies 
increase with decreasing temperature. RPA also predicts 
a fictitious negative thermal expansion coefficient at this 
density as shown in Fig. @] (b). In high densities, how- 
ever, agreement of RPA with simulation is excellent for 
all temperatures down to the freezing temperatures as 
shown in Fig. [4] (c)-(f). This is very surprising because 
the temperatures in these figures are far below Tth > where 
RPA fails to describe overall shapes of S(k) as discussed 
in the previous subsection. 

We also calculated T a< spA by solving <9Prpa/<9T = 
and plotted in Fig. [TJ (dotted line). Agreement of T a RPA 
with the simulation result is perfect except for the vicin- 
ity of the re-entrant melting region. The asymptotic ex- 
pression of T^ppa which is valid at high densities can be 
written as (see Appendix) 



te^aHPA 3 / 2 3 

; = 7T ' per exp 



-pa' 



2/5 



.(17) 



This asymptotic expression works well down to p* <~ 1 
as in the case of T t h- Interestingly, the density exponent 
2/5 in this expression is smaller than 2/3 for T t h and 
T m (see Eqs. © and (JTHJ). Due to this difference, T a 
monotonically deviates from T m as the density increases 
and it eventually becomes larger than T t h in the high 
density limit (not shown). 

We made similar comparison for the potential energy 
in Fig. [7] Fig. [7] (a) and (b) show the temperature de- 
pendence of the potential energy in the fluid phase (filed 
circles) and crystalline phase (open circles) at two den- 
sities. The uniform part of u is temperature indepen- 
dent (see Eq. ©), although it dominates the net values 
of u. In order to see the temperature dependence of u 
more clearly, the fluctuation part Am are shown in Fig.[7J 
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FIG. 7: The temperature dependence of u and Alt at p* — 
0.33 ((a) and (c)) and p* = 2.01 ((b) and (d)). Filled and 
open circles represent the simulation results for the fluid and 
crystalline (bec) phase, respectively. Solid lines are the RPA 
results for the fluid phase. 
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FIG. 8: Temperature dependence of the ratio of RPA to sim- 
ulation values of the fluctuation part of the potential energy 
(left) and the pressure (right) for p* = 0.33 (O), 1-00 (□), and 
2.01 (A). The range of temperatures in these figures is much 
lower than the corresponding Tth- (T^/T^ are 1.48 x 10 -2 , 
3.14 x 10" 3 , and 2.55 x 10~ 4 for p* = 0.33, 1.00, and 2.01, 
respectively.) 



(c) and (d). One observes that agreement of the simu- 
lation results with RPA is far better at high densities; 
At p* = 0.33, RPA underestimates Am and discrepancy 
from the simulation data are larger than the energy gap 
between the fluid and crystalline phases. At the higher 
density p* — 2.01, the discrepancy is less 10% even at 
the freezing temperature. 

In order to quantify the accuracy of RPA for both the 
potential energy and pressure, we plot the ratios of Am 
and AP of simulation to those of RPA in Fig. [5] against 
the inverse temperatures for several densities. The ratios 
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tance decouples with the length scales which dominate 
thermodynamic quantities of GCM. This is the reason 
why RPA remains the excellent approximation to predict 
thermodynamic quantities even far below the threshold 
temperatures. 
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FIG. 9: The integrand of Eq. (O (shaded areas), the inte- 
grand where S(k) is replaced with SrpaCc) (solid lines) and 
S(k) (broken lines) for (a) T* = 7.9 x 10" 3 and p* = 0.33, 
(b) T* = 5.0 x 10~ 4 and p* = 1.00, (c) T* = 3.2 x 10~ 5 and 
p* = 1.65 and (d) T* = 7.9 x 10~ 6 and p* = 2.01. 



decrease rapidly as the temperature decreases at the low 
density p* = 0.33, whereas, in the high density p* = 2.01, 
the ratios for both the potential energy and pressure re- 
main to be more than 80% for the whole range of tem- 
peratures down to the melting temperature. 

Given that RPA poorly describes even the qualitative 
behaviors of S(k) and g(r) at the high density and low 
temperature regimes, it is very surprising and counterin- 
tuitive that RPA is excellent at predicting quantitatively 
thermodynamic quantities u and P. In order to ratio- 
nalize this puzzling facts, we look again the integral ex- 
pressions of the thermodynamic functions, Eq. ([9]) and 
(fT2"]) . These expressions show that both Am and AP are 
expressed in terms of the integral of S(k) multiplied with 
the pair potential over the wavevectors. In order to see 
which length scales dominate the integrand, we show the 
integrand of Am in Eq. @ with simulated S(k) in Fig. [5] 
(shaded area), along with those obtained using 5rpa(&) 
(solid line). At the low density p* — 0.33, the integrand 
with simulated S(k) is peaked at the peak position of 
S(k) (broken line). This implies that Am at this den- 
sity is dominated by the contribution at the interparticle 
distance, just like ordinary atomic fluids. The integrand 
obtained using RPA fails to account for this peak struc- 
ture (solid line). However, with increasing density, the 
peak position of the integrand shifts to smaller wavevec- 
tors than the first peak position of S(k). Concomitantly, 
agreement of the integrand obtained from simulation and 
RPA becomes better and better. This agreement origi- 
nates from the fact that RPA can account excellently for 
the low wavevector behavior of the static structure factor 
as shown in Fig. [SJ At very high densities, the particles 
start overlapping and the characteristic interparticle dis- 



In this paper, we have presented a detailed analy- 
sis of thermodynamic and structural properties of the 
high density one-component GCM. Special emphasis has 
been put for static properties of the fluid phase. First, 
the solid-fluid phase boundary of the system is care- 
fully evaluated up to the unprecedentedly high density 
p* = 2.4. Our result confirmed the scaling conjectured 
by Stillinger for the freezing and melting temperatures, 
log!/, logT m cx -p 2/3 , at p* > 1.2. The potential en- 
ergy difference between the crystalline and fluid phases 
was shown to be linear in the freezing temperature and 
the entropy difference is almost constant at high densi- 
ties, which verifies the assumption which Stillinger's ar- 
gument is based upon. The thermodynamic and struc- 
tural properties of GCM in the fluid phase are analyzed 
in detail for a wide range of temperature and density. 
The potential energy u, the equation of state P, the 
static structure factor S(k), and the radial distribution 
function g(r) were evaluated by simulation. We com- 
pare the simulation results with the RPA results. In the 
high temperature regime, RPA provides almost perfect 
description for both thermodynamic quantities and the 
structural factor S(k). RPA is rather poor at predicting 
g(r) at r ~ 0. Threshold temperature T t h below which 
RPA fails to describe g(r = 0) is relatively high. In 
the high density and low temperature regime, RPA fails 
to capture the peak structure of S(k) even qualitatively, 
whereas it predicts correctly the low fc's behavior up to 
just below the wavevector at which the first peaks are 
located. Despite of poor performance of RPA at describ- 
ing the structural properties, RPA successfully describes 
thermodynamic quantities such as the potential energy 
and pressure at high densities. Agreement of RPA with 
simulation results systematically improves as the density 
increases even near the phase boundary. The temper- 
ature below which the thermal expansion coefficient be- 
come negative is also accurately calculated from RPA. By 
scrutinizing the role of the microscopic structure of par- 
ticles in the potential energy and pressure, we concluded 
that the surprising success of RPA is originated from the 
decoupling of the length scales which dictate the thermo- 
dynamic quantities and the interparticle distance. This 
decoupling is attributed to the mild and long-ranged re- 
pulsive tails of the pair potential of GCM. The fact that 
RPA is an excellent approximation even at the vicinity of 
the phase boundary, or even at the supercooled regime, 
at high densities, hints that the mean-field description is 
valid for the high density GCM and may play a crucial 
role to understand (glassy) dynamics let alone thermo- 
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FIG. 10: The validity of the asymptotic expressions for Ty, 
and T a rpa- Open circles and solid line are RPA values and 
its asymptotic expression, Eq. (|16[) . of Ti^, respectively. Open 
squares, dotted line, and dashed line are RPA values and its 
two asymptotic expressions, Eqs. (|A4|) and (|17[). of T^rpa, 
respectively. 



dynamic properties [35j, [36[ . 
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Appendix A: Derivation of Eqs. (11G[) and (|1T[) 

In this appendix, we derive the asymptotic expressions 
of T t h and Tq.rpa at high densities b y u sing the asymp- 
totic expansion of the polylogarithm |42| 

LU-z) = -^^y + 0((-\o gx y- 2 ), (Al) 

where T(x) is the gamma function. 



At T = T th , Eqs. (S) and (JM) lead to 
^* - ^2 (1°8^) 3/2 + ^ ((logx- t/l )- 1/2 ) 



(A2) 



where = TT 3 ^ 2 p* /T^ h . When p* is sufficiently large, we 
can neglect the third term on the left hand side, leading 
to the asymptotic expression of T t h of Eq. (|16[) . Fig. QJJ] 
shows that this asymptotic expression is very accurate 
down to p* ~ 1. Likewise, at T = T Qi rpa, Eqs. (fl"4"|) and 
(|AT|) lead to 



— ^ (logXR PA ) 5/2 + 77~"~2 (logXRPA 
157T^ ^ Z 



4 
37 2 



.3/2 



(A3) 

+ o((lo ga; RPA) 1/2 ) =0, 
where zrpa = n 3 / 2 p*/T* RP a- If wc neglect the third 
and fourth terms on the left hand side, we obtain 



RPA 



= tt 3 / 2 



p exp 



15^V x2/5 



(A4) 



However, Fig. [TU] shows that this expression is not accu- 
rate even at very high densities. For sufficient accuracy, 
we have to keep the third term on the left hand side of 
Eq. (|A3|) . If we neglect only the forth term on the left 
hand side, Eq. (|A3[) becomes 



,.- 2 /BjL fl 3/ 2 _ 1 = 0) (A5) 



15tt 2< 



a 5/2 



where s = p*~ 2 / 5 log(7r 3 / 2 p*/T*). Since the second term 
on the left hand side of this equation is small due to 
the factor p* -2 / 5 , we can expand the solution as s = 

50 + p*~ 2 / 5 si + • • ■ and can solve the equation in each 
order of density. The first order solution including sq and 

51 is 



15tt 2 \2/5 



2p* 



-2/5 



(A6) 



which leads to the asymptotic expression of T Qj rpa of 
Eq. p7|) . Fig. [TU1 shows that this expression is very ac- 
curate down to p* ~ 1. 
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